• Splines

• Regression

• Adaptive

• Multivariate

• They are flexible lines/curves that are combined together

• Individual splines components are generally relevant to only part of the space/range of a feature or a set of features

• They have continuity constraints (joined together - no jumps), and smoothness constraints (e.g. up to second derivative continuity at a join)

library(tidyverse)

require(ISLR)
?Wage #for more details on the dataset

# Exclude very high salaries
Wage <- Wage %>% filter(wage < 250)
p <- ggplot(Wage, aes(x=age,y=wage)) + geom_point(colour='blue')
p

• MARS tries to use as small number of basis functions as ‘nessesary’
• It is adaptive to the data, using only the features and flexibility required
• Uses a build phase (adding basis functions) and a pruning phase (removing basis functions)
• This applies to the features and the knots (breakpoints) used within a feature

• Splines can be fit with the same machinery as multivariate linear models, these use a ‘design’ or ‘model’ matrix

• Each column of your model matrix corresponds to a basis function, more on this later

• Interactions are modelled using tensor product. This is equivalent to including a column with element-wise multiplication of other columns (same as in standard multivariate linear regression)

• To fit a linear piecewise model we can use the machinery of linear models

• We use spline basis functions

• Run a feature through the basis function. This is the same idea as when you would square a feature, or log it, truncate to max value before using in linear model

• Use the resulting features as the columns of the model matrix \(\bf{X}\) in your linear model ‘\(\bf{y = a + Xb + \epsilon}\)’

• Class attributes can be added to model matrix with usual one-hot encoding

• It is a generalisation of linear modelling, MARS assumes Gaussian errors

• The type of basis functions used in MARS are of the form

• \(\bf{(x-t_{x3})_+}\)

• \(\bf{(t_{x3}-x)_+}\)

and interactions such as:

• \(\bf{x \cdot y}\)
  
• \(\bf{x\cdot(y-t_{x4})_+}\)
• \(\bf{(x-t_{x3})_+\cdot(y-t_{y2})_+}\)

where \(\bf{x}\), \(\bf{y}\) are features and \(\bf{t_{x3}}\) etc. are knots (actual values in the range of the feature).

• Start with a base model (intercept)

• Scan through all candidate basis functions (all features and all knots)

• At each stage we consider the product of basis functions in the model with candiate basis functions

• One restriction is that a feature can appear only once in a term (interactions ok, higher order powers are not)

• Usually restrictions are put on the degree of interactions

• Add the basis functions pair which most reduces the loss (RSS)

• Continue until a stopping criterion is TRUE (e.g. max number of basis functions, reduction in RSS is very small, numerical stability)

• Remove basis functions one-by-one, each time choosing the one with the least reduction in Residual Error (the most useless)

• The ‘Best’ model comapring all reduction phase models is found by (default) GCV Generalised Cross Validation which is similar to AIC or BIC. GCV prefers lower RSS but penalises more complex models over less

• GCV stands for Generalized Cross Validation, a perhaps misleading term because no cross-validation is actually performed

• One can also use a test set to pick the best model, or real cross validation - depending on implementation

Important point

• The term “MARS” is trademarked and licensed exclusively to Salford Systems: https://www.salford-systems.com. We can use MARS as an abbreviation; however, it cannot be used for competing software solutions.

• The only official MARS is from Salford Systems

• This was published as a package in Splus before the trademark restrictions came in

• Coded in C and S

• Slightly different algorithm to real MARS

• Allowed for multiple responses (polychotomous), hence classification

• Although open-source, Salford System licenced the commerical use of the code for 10 years

• The package was ported in R, eventually the functionality was incorporated into a more general spline package called polspline

• MDA is a package that originally implemented MARS in R, authors Hastie & Tibshirani

• The package is focussed on Mulivariate Discriminant Analysis but has full MARS functionality

• Earth package (in R and Python) is probably the best version to use

• Implements MARS

• Does Polycotomous Regression (so Classification)

• Implements many Loss functions using the machinery of GLMs

  • Binomial
  • Poisson
  • Gamma
  • Tweedie
  • inverse.gaussian

This means successes; counts; insurances losses; and randomly large positive amounts, etc. can be modelled

library('tidyverse')
library('earth')
library('caret')

input_dir <- '/Users/Martin/Google Drive/Projects/MARS'
data_file <- "housing.data"

housing <- read_table(file.path(input_dir, data_file), col_names = FALSE) 
colnames(housing) <- c('CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT', 'MEDV')

set.seed(92835)
training_idx <- createDataPartition(housing$MEDV, p=0.7, list=FALSE)
train_data <- housing %>% slice(training_idx)

emars.fit <- earth(MEDV ~ . , degree=2, trace= 3 , data = train_data)

• Benchmark ML dataset

• 28 x 28 grid of pixel intensity 0 - 255

• Training set: 60,000 digits

• Test set: 10,000 digits

• http://yann.lecun.com/exdb/mnist/

• Degree 1, multiclass Earth; 10 of 784 predictors, 21 terms: 66%

• Degree 3, multiclass Earth; 20 of 784 predictors, 41 terms :74%

• Degree 4, multiclass Earth; 34 of 784 preds, 77 terms: 78%

• Degree 5, multiclass Earth; 14 of 784 predictors; 27 terms : 64%

• Preprocessed the data with de-skewing

• Used 10 model ensemble of one-versus-rest binomial loss EARTH models

• Balanced data in each single model through down-sampling
  

• Maximum degree of interaction permitted: 5

• Test set accuracy of 95.4% - not bad

• But don’t use a model like this for image analysis

• https://archive.ics.uci.edu/ml/datasets/Wine

• These data are the results of a chemical analysis of wines grown in the same region in Italy but derived from three different cultivars. The analysis determined the quantities of 13 constituents found in each of the three types of wines.

1 Alcohol 2 Malic acid 3 Ash 4 Alcalinity of ash 5 Magnesium 6 Total phenols 7 Flavanoids 8 Nonflavanoid phenols 9 Proanthocyanins 10 Color intensity 11 Hue 12 OD280/OD315 of diluted wines 13 Proline

• Number of Instances - Class 1 has 59; Class 2 has 71; Class 3 has 48

• Train/test : 70% / 30%

library('earth')

...        

emars_fit <- earth(x = train_feat, y = train_label_one_hot, degree=2, trace= 4 )
e_test_pred <- apply(predict(emars_fit,newdata = test_feat),1,which.max)
e_confusion_tab <- table(e_test_pred,test_label)

library('polspline')

...

pmars.fit <- polymars(resp=train_label,pred = train_feat, verbose=TRUE, classify = T)
p_test_pred <- predict.polymars(pmars.fit, x = test_feat, classify = TRUE)
p_confusion_tab <- table(p_test_pred,test_label)

• MARS by way of EARTH is a handy package for whipping up a linear/GLM model

• Anywhere you use a linear model, MARS can quickly provide

  • Feature selection
  • Feature Engineering
  • Interaction selection
  • and a decent model

• Many extra features we didn’t cover like fine grained model search such as: terms not allowed to interact, start model etc.

• Good documentation

• Works with mixed data-types, continuous and class (one-hot encoding)

• Simple to set up, not many configuration options

• Works across a number of GLM types (Tweedie, Gamma, Binomial, …) as well as Gaussian (MSE)

• Can get you most, if not all of the way there, and without much fuss

• or it can provide a good initial benchmark model if you are going to invest a lot of time in other forms of models